arXiv:math/0512347vl [math.NA] 14 Dec 2005 


Single Exponential Approximation of Fourier 

Transforms 

Patrick McLean* 

(31 November 2005) 


Abstract 

This article is concerned with a new method for the approximate 
evaluation of Fourier sine and cosine transforms. We develope and 
analyse a new quadrature rule for Fourier sine and cosine transforms 
involving transforming the integral to one over the entire real line and 
then using the trapezoidal rule in order to approximate the trans¬ 
formed integral. This method follows on from the work of Ooura and 
Mori, see [Hj and [Ij 

A complete error analysis is made using contour integration. An 
example is examined in detail and the error is analysed using residues 
and the saddle point method. The method we have developed is char¬ 
acterised by its simplicity and single exponential convergence. 
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1 Introduction 

In this paper we consider the numerical approximation of Fourier sine and 
cosine transforms, that is, integrals of the form 


fc{t) = 

poo 

/ f{x) cos{tx)dx, 

Jo 

( 1 ) 

fs{t) = 

poo 

/ f{x) sm{tx)dx. 

Jo 

( 2 ) 


There exist extensive tables of Fourier transforms, see, for example, j2]. 
Nonetheless there is need for the numerical approximation of Fourier trans¬ 
forms. We focus on fs{t) and present a new method involving transforming 
the integral m to one over (—exo, exo) and then using the trapezoidal rule. 

In section 2 we review the trapezoidal rule on (— cxd, cx))and sources of 
error in its use. Quadrature rules over intervals other than (—exo, cx)) arise 
from the use of transformations. The double exponential and sine methods 
of quadrature arise in this way, see PI and [SI, respectively. 

In section 3 we introduce quadrature methods for Fourier sine integrals 
based on transformations of a certain sort. Double exponential versions were 
introduced in jS] and [7|. Here we present a single exponential version to¬ 
gether with useful asymptotic estimates of its behaviour. We also note that 
the midpoint rule may be used for Fourier cosine transforms. 

In section 4 we present an example and analyse the discretisation error 
using residues and the saddle point method. Here we see that the rate of 
convergence of our method is determined by the proximity of singularities 
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of the integrand to the interval of integration (0, cxo). This phenomena is 
common in quadrature, see jS], for example . 

In section 5 we look at the truncation error and compare our method with 
the double exponential method of Ooura and Mori. 


2 The Trapezoidal Rule and its Error 

Given a function F{u) dehned on (— cxo, cx)) with integral 

/ OO 

F{u) du, (3) 

•CO 

and a real number h > 0 we dehne the trapezoidal approximation to I by 

CO 

n = hJ2 Fikh), (4) 

k = — OQ 

This rule is the basis of Sine methods extensively developed by Stenger and 
can be obtained by integrating a Sine function interpolant to F{u), see j2]. 


2.1 Discretisation Error of the Trapezoidal Rule 

The quantity I — Th is referred to as discretisation error. We shall use the 
error representation of Donaldson and Elliott jS]. If we set A = 0, a = 0 and 
u = 1/h then equations (6.1) and (6.4) of jH] read 


'h/i(w) 

*hh(w) 


) 7iexp{ifw), > 0 

TT exp(—i|tc), < 0 

-sm{-w), 


(5) 

( 6 ) 


respectively. Thus, from jHl eqn 2.4] we have the error representation 


I-n = 


1 f ^fe(w) 
27r* Jc ^h{w) 


F{w) dw, 


(7) 


where G is a positively described contour enclosing the zeroes of <h/i(w)and 
avoiding any singularities of F{w). 
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2.2 Truncation Error of the Trapezoidal Rule 

In practise, the (infinite) trapezoidal snm Th is trnncated at some ±n to give 
the finite trapezoidal sum Tn^h- 


n 



( 8 ) 


k=—n 


as an approximation to /. This introduces an error Th — Tn^h that we refer 
to as truncation error. 

Thus, the overall error can be seen as comprising two sources: 


(9) 


I — Tn^h — I — Th + Th— Tn^h 


It is essential that the discretisation error and the truncation error match in 
order to achieve a suitable rate of convergence. 

3 Trapezoidal rule methods for Fourier sine 
transforms 

In this section we propose a method of approximation of Fourier integrals 
which is highly accurate and applicable to a wide range of functions includ¬ 
ing ones with singularities and slow decay. Our approximation to fs{t) is 
obtained by making a particular change of variable from (0, cxo) to (—cxo, cxo), 
followed by an application of the trapezoidal rule. For fdt) we apply the 
same change of variable, but use the midpoint rule. 

In order to approximate fsit) we take the integral 0 and introduce the 
change of variable 



—oo < u < oo, 


( 10 ) 


where m > 0 and 0 : (—cx),cx)) —> (0,cx)) is a function satisfying 


0(m) ~ 0 as M —>• —oo, 

(f){u) ~ M as M —> -l-CXD. 


( 11 ) 

( 12 ) 


Thus, fs{t) can be represented by the integral 



( 13 ) 
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where the function F^iu) is dehned by 

TD 771 

FmH = /(y0(n))sin(m0(n))y0'(n). (14) 

We now apply the trapezoidal rule with stepsize h = ^ to (HI to give the 
approximation 

OO 

Tm = - y2 Frnik-). (15) 

m m 

k=—oo 

The asymptotic behaviour of 0 (m) at ±cxd allows us to truncate this series 
without incurring too large an error. We shall provide further details on this 
in the last section of this paper. 


3.1 Midpoint rule for Fourier cosine transforms 

For completeness we note that fc{t) may be efficiently evaluated by making 
the transformation x = ^4>{u) and then applying the midpoint rule 

OO ^ 

M, + -)A), (16) 

k=—oo 

with stepsize h = ^. The analogous error representation is 

I - M, = ^ dw, (17) 

21X1 Jc ^h{w) 

where and ^h{w) are given by 

f-z7rexp(ifM;), > 0 

IZ7rexp(—z^tn), < 0 

cos{^w), (19) 

respectively. Further details including an example are provided in M- 


Tft(w) = 
4>ft(w) = 


3.2 Double exponential transformation I 

In their hrst paper [HI eq 11] Ooura and Mori introduce the transformation 


01 (u) 


1—exp(—X sinh u) ’ 

1 


ti ^ 0 

M = 0 


( 20 ) 


for K > 0. They make the choice of iF = 2ti. 
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3.3 Double exponential transformation II 


In their second paper [3 eqn 3.3] Ooura and Mori introduce the transforma¬ 
tion 


02(m) = 


U 


1 — exp(—2 m — a{l — e““) — /3(e“ — 1)) 
where the parameters are given by 


( 21 ) 


a = /3/\/l + Af log(l + M)/(2ir), /3=j' 


( 22 ) 


3.4 Novel Single Exponential Transformation 

We dehne the transformation 0 : (—oo, cx)) —t (0, cx)) by 

0(m) = log(e“ -|- 1), —cx < u < oo, (23) 


with derivative 


Note that as a function of the complex variable w 
points at w = {2k + l)i7i for an integer k. 


(24) 


u + iv, (l>{w) has branch 


3.4.1 Asymptotic behaviour of transformation 

We shall need to know the asymptotic behaviour of 0 (m) as m —t ±cx. The 
Taylor series of log(l -|- x) around a; = 0 is 

log(l -1- x) = a; -I- O(x^), for |x| < 1, (25) 

see P eqn 4.1.24]. Thus, we have that 

0(m) = 0-Fe“ + O(e2“), (26) 

as M —t —oo, and 

0(m) = u + log(l -1- e““) =u + e““ -h 0(e“^“), (27) 

as M —t X. 
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3.4.2 Asymptotic behaviour of inverse 

The inverse 0“^ : (0, cx)) —>■ (—oo, cx)) is defined by 

= log(e'^ ~ l)^ 0 < x < cx). (28) 

We shall determine the asymptotic behaviour of w = (j)~^{z) as 2 ; —> 0. 
From (PH|) we have 

- 1 

w = log(e^ — 1) = log^; + log-. (29) 

‘Z/ 

On expanding the last term about ^ = 0, we arrive at 

0“^(z) ~logz + | + 0(z^), (30) 

as z —> 0 , argz 7 ^ tt. 


4 A Certain Integral 


In this section we shall consider the integral 


1 



sin tx ^ 
{x - a )2 + 62 


(31) 


for t > 0, — CX) < a < 00 and 6 > 0. We shall estimate the discretisation error 
using residues and the saddle point method. We shall see that the discreti¬ 
sation error depends on the proximity of the singularities of the integrand 
a ± to the interval of integration ( 0 , cx)). 


4.1 Evaluation in Terms of Trigonometric Integrals 

We follow the standard reference work P §5.2] and define the trigonometric 
integrals as follows. For ^ G C define the sine integral Si(z) as 

Si(z) = / - dt. (32) 

Jo ^ 

The sine integral Si is an entire function. A commonly used notation is 
si(z) = f - Si(z). 

For z E C such that | argz| < tt define the cosine integral Ci(z) as 

POS f — 1 

Ci(^) =7 + logz+ / ^ - dt. (33) 

Jo i 
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The cosine integral Ci has a branch cut discontinuity along the negative real 
axis. 

The functions Si and Ci occur as Fourier sine and cosine transforms. 
Specifically, for | arga| < vr and y > 0 we have that 



cos{xy) 
a + X 
sin (a:?/) 
a + X 


dx 

dx 


— si{ay) sin(a?/) — Ci{ay) cos{ay) 
Ci{ay) sm{ay) — si{ay) cos{ay), 


(34) 

(35) 


see 0 §1.1 (9)] and 0 §1.2 (10)], respectively. 

Using partial fractions together with (ld4D and (jnsi) it is possible to show 
that 


bl = sin at sinh bt^ Ci{—a + ib)t + cos at cosh a + ib)t (36) 

— cos at sinh bt^ si{—a + ib)t — sin at cosh bf^ Si(—a + ib)t). (37) 


4.2 Discretisation error 

In order to implement our method we make the change of variable x = y0('u) 
and apply the trapezoidal rule with stepsize — to give the approximation 


T = 

m 



where the function Fm{u) is given by 


sin(m0(M)) Y(;/)'(k) 

^ (f0(M) -a)2 + 62- 


(38) 


(39) 


We shall denote the pole of F^iw) in the upperhalf plane closest to the real 
axis by wq and we note that wq is also a pole of Fm{w). 

Now, from (ED we have that 


= ^ dw, (40) 

2m Jc ^m[w) 

where ^/^(tc) and $m('ia) are given by® and respectively, and where C 
is a positively described closed contour going between the real axis and the 
poles wq and wq as depicted in hgure [U 







Figure 1: The contours C and C 

If we deform the contour to contain the poles Wq, FJo we have by Cauchy’s 
theorem 

I -Tra = -Res{^Fm-,wo) -Res{^Fm]w^) + ^ dw, 

2tII Jc> ^m{w) 

m 

where C' is a positively described contour going between the poles wq, wq 
and the singular points of (j){w) as depicted in £gure[TJ We shall denote the 
two residue terms by Rm and the integral term by S^- 


4.3 Contribution to the error from Poles 

In this section we consider the evaluate the residue terms in (EH)- 

Now, on using L’Hopital’s rule, the residue of Fm{w) at wo is given by 

^ sin(m0(ii;)) Y0'(w) 


Res{Fm',wo) = lim (ic — tco) 

w— 

sin(m0(tco)) 


w^wo [ ^(f){w) — aY + IF 


2 (f0(n;o)-a) 


Now, we use the fact that m(j){wo) = (a + ib)t to give 

sin((a + ib)t) 


Res{Fm-,wo) = 


2ib 


Similarly, we have that the residue of Fm{w) at Wq is given by 

sin((a — ib)t) 


Res{F^;wo) = 


-2ib 


(42) 

(43) 

(44) 

(45) 
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Combining these residues with the dehnitions of and ^rn{w) we 

have that 


R 


m 


71 exp{imwo) sin(a + ib)t ti exp{—imWo) sin(a + ib)t 
— sm{imwo) 2ib — sin(zmwo) —2ib 

71 sin(a + ib)t ^ sin(a — ib)t 
b \_1 — exp{—2miwo) 1 — exp{2miwo) 


(46) 

(47) 


We now substitute wq = uo + ivo and wq = uq — ivo and rearrange to give 


^ 71 [exp(—2muo) — cos(2mMo)] sin(at) cosh(6t) + sin(2mMo) cos(at) sinh(6t) 

b cosh(2muo) — cos(2mMo) 

(48) 

Now, on neglecting the cosine term in the denominator and the exponential 
terms exp (—27?wo) in the numerator we have that 


27r — cos(2mMo) sin(at) cosh(6t) + sin(2muo) cos(at) sinh(6t) 
b exp(2muo) + exp(—2muo) 


(49) 


On neglecting the term exp (—2muo) in the denominator, we obtain 
27r 

Rm ~ -g- exp(—2muo) [— cos(2mMo) sin(at) cosh(6t) + sin(2mMo) cos(at) sinh(fet)] 

(50) 

We shall only be interested in the asymptotic behaviour of tco as m —>■ cx), 
which we can obtain from dsni); 


, {a + ib)t ia + ib)t _ 2 , 

Wq = log ^ ^ ^ ^ ^ + 0{m 2), 


m 


2m 


(51) 


as m —^ oo. Now since log(re*®) = log(r) + iO for r > 0, 0 < 0 < tt, it follows 
that the imaginary part vq of wq has asymptotic behaviour 


farctan(^) + |h + o(^), a > 0, 
no = h + ^, a = 0 (52) 

U-arctan(^) + ^ + 0(;^), a < 0. 

as m —>• cx). 

Thus using the asymptotic estimates dSDD and (IS2D we have the estimate 


[ C exp (—2m arctan(^)) , 

a > 0, 


-Rm ~ < Cexp (-2m|)) , 

a = 0 

(53) 

exp (—2m(7r — arctan(^))) , 

a < 0, 



as m —> cx), where C represents a number independent of m. 
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4.4 Contribution to the error from the saddle points 

In order to estimate the contribution to the error I — from the integral 
Sm in m we write its integrand as exp{p{w))q{w) where 

\\j ('}n\ 

p{w) = log(—^ sin(m(;/)(w)) (54) 

Tfl 771 

q{w) = (55) 

From the dehnitions of '^rn{w) and <hm(w) we have that 

p{w) = log (—7r(cot(mtc) + i) sin(m0(tc))), '^{w) > 0. (56) 

The saddle points of p{w) are solutions of the equation p'{w) = 0. We shall 
be interested in the one in the upper half plane closest the real axis and we 
shall denote it by Wi. The point uq in the lower half plane is also a saddle 
point of p{w). 

To hnd Wi we must solve 

p'{w) = —m{cot{mw) — i) + cot {m(j) (w)) mcl)' (w) = 0, (57) 


that is, 

(j)'{w) cot(m0(tc)) = cot(mtc) — i. (58) 

Solving this equation is problematic. Possible strategies are i) obtaining an 
asymptotic estimate for wi in m similar to the asymptotic estimate tco given 
by dSII), or ii) given values of m numerically solving for Wi. Lacking the 
former we use the latter. 

From m eqn2.7.2] we have 


d'm(w) 


27r 


9 • / ff) t \^rn{w)dwr^—J (59) 

2m Jc+ <hm(w) 2,m y p"{wi) 


where is the upper half of C' and where 

« = exp(y - ^ arg(/(u;i))). 


(60) 


The contribution to the error from the saddle point at Wi is the conjugate of 
the estimate in Hence, we have the estimate 


~ 23? 


27r 


-ae 


p(wi) 


27ri y p''{wi) 
as an estimate of the integral over C' (see (HOD) 


q{wi) 


(61) 
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a = —1 

a = 0 

a = 1 

m 

I-Tm 

Rm 

Sm 

I-Tm 

Rm 

Sm 

1 

'— 1 

Rm 

Sm 

1 

8.60E-3 

2.20E-2 

-1.40E-2 

-2.50E-2 

O.OOEO 

-1.40E-2 

-2.50E-2 

O.OOEO 

-1.40E-2 

2 

-1.60E-4 

-2.60E-4 

-2.80E-4 

-1.60E-3 

-1.80E-3 

-5.60E-4 

-1.60E-3 

-1.80E-3 

-5.60E-4 

3 

7.90E-6 

2.30E-6 

-4.80E-9 

-6.70E-5 

-6.70E-5 

-4.10E-9 

-6.70E-5 

-6.70E-5 

-4.10E-9 

4 

-3.40E-8 

-2.00E-8 

-7.00E-9 

9.40E-6 

9.40E-6 

-2.50E-8 

9.40E-6 

9.40E-6 

-2.50E-8 

5 

-4.90E-9 

l.OOE-11 

7.40E-9 

1.60E-7 

1.50E-7 

7.80E-9 

1.60E-7 

1.50E-7 

7.80E-9 

6 

-4.50E-11 

1.80E-12 

1.70E-12 

-8.20E-9 

-8.30E-9 

1.90E-12 

-8.20E-9 

-8.30E-9 

1.90E-12 

7 

4.90E-12 

4.40E-15 

9.40E-12 

-6.50E-10 

-6.60E-10 

9.80E-12 

-6.50E-10 

-6.60E-10 

9.80E-12 

8 

1.20E-13 

-8.70E-17 

-4.60E-14 

-3.10E-11 

-3.20E-11 

-3.60E-14 

-3.10E-11 

-3.20E-11 

-3.60E-14 

9 

-4.30E-15 

-1.20E-18 

-1.20E-14 

-1.40E-12 

-1.40E-12 

-1.20E-14 

-1.40E-12 

-1.40E-12 

-1.20E-14 

10 

-l.lOE-16 

-l.lOE-20 

1.50E-16 

-5.30E-14 

-5.40E-14 

1.40E-16 

-5.30E-14 

-5.40E-14 

1.40E-16 


Table 1: The error / — Tm and contribution from residues Rm and saddle 
points Sm for a = —1, 0,1 


4.5 Discretisation error 


In this section we present results comparing the discretisation error / — Tm 
with the contribution from the poles Rm and the contribution from the saddle 
points Sm- 

To evaluate / we use equation (EZD and to evaluate Tm to within machine 
precision we sum the series for sufficiently high values of n. Typically taking 
n = AwS. Furthermore, we use the representation (P|) for Rm and the 
estimate (jnH) for Sm- 

We evaluate the three quantities / — Tm, Rm and Sm for m = 1,4, ...10, 
a = —1,0,1, t = 1 and 6=1 and presenting the results in Table [Hand we 
plot their base 10 logarithms in Figure |21 

We observe that for a > 0 the error I — Tm is determined by Rm and Sm 
is negligible in comparison. Furthermore, we observe that 


I — Tm ^ Rm ~ 


Cexp (—2m arctan(^)) , 
Cexp (-2m|)) , 


a > 0, 
a = 0. 


( 62 ) 


While for a < 0 the error J—Tm is determined by Sm while Rm is negligible 
in comparison. Furthermore, we observe that 


I - Tm ^ Sm ^ C'exp (-vrm)). 


(63) 
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5 Truncation error and comparison with dou¬ 
ble exponential case 


5.1 Matching Truncation and Discretisation Errors 

We now estimate the error introduced by truncating the inhnite series Tm to 
the hnite series dehned by: 


T = 

n,m 



(64) 


This error is bound by 

\Tm-Tn,m\ < “ \ f {m(j){k—)) sm{m(j){k—))m(j)'{k—)\ (65) 

m m mm 

\k\>n 


We we have that \f{x) \ is bounded by some constant Cf. Also, we have that 


, TT e ^ 

^\k-) = < 1. 

m e rn + 1 


( 66 ) 


Thus, we must investigate the asymptotic behaviour of sin(m0(fc^)) as k 
±cx). First, using we have, as k ^ —oo, that 


vr 


sin(m0(A;—)) = sin(me'‘™ + )) 

m 

= + 0(e^^^). 


By we have, as k —> oo, that 


TT 


sm(m(j)(k —)) = sin(/c 7 r + me ^ + 0(e ”*)) 

m 

= (—1)^ sin(me“^^ + 0(e“^^^)) 
= (—l)^me“^^ + 0(e“^^™). 

Thus, from (ESI) we have that 

OO 

\T^-Tn,^\ < 2- V C'y(me-'^^+0(e-2"^)) 

k 

< 27iCf- 


k=n+l 

S ^ ^ m 


1 — e 


< 27tC 


r 


(67) 

( 68 ) 

(69) 

(70) 

(71) 


(72) 

(73) 

(74) 
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Figure 3: Our method vs Ooura and Mori’s First method for Ii 


Thus, if the discretisation error is, say, 

/ - ~ (75) 

then we equate the exponents in the previous two equations to determine the 
dependence of m and n, that is. 



This choice of m results in an overall error 

I Tyi^jYL ^ Cy,/^C 


(76) 


(77) 


5.2 Single exponential vs Double exponential 

In Figure El we present the errors in approximating the integral 


Ii = 


sm X 


X 


dx 


(78) 


using the hrst method of Ooura and Mori and using our method. Clearly 
our method results in slower convergence than theirs. 
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